namespace hms::swe
{

template<typename MeshType>
template<BlockLocation loc>
scalar SolverThread<MeshType>::timestep(
	const Block<loc>& block,
	const Ref<const ArrayXXs>& dBlock,
	const Ref<const ArrayXXs>& vBlock
){
	const MeshType& mesh { solver().mesh() };

	auto courantFactor = [&](
		[[maybe_unused]] Index y
	){
		if constexpr ( std::is_same_v<MeshType, UniMesh> ){
			return 4 * mesh.cells().area() / mesh.cells().circumf();
		} else if constexpr ( std::is_same_v<MeshType, RectMesh> ){

			/* if areas or circumferences are precalculated,
			 * then mesh.cells().areas() and mesh.cells().circumfs()
			 * reference dense arrays rather than expressions. Then, each
			 * block row can be mapped. */
			// #define MAP(ARG) Array1XsMap{ &(ARG[0]), 1, block.nx() }

			// auto getReadRow = [&]( const auto& arr ){
			// 	return block.template writeBlock<1>(arr).col(meshRow);
			// };

			// #undef MAP_IF_DENSE

			// const auto& areas {
			// 	#ifdef RECTMESH_PRECALC_AREAS
			// 	MAP( getReadRow( mesh.cells().areas() ) )
			// 	#else
			// 	getReadRow( mesh.cells().areas() )
			// 	#endif
			// };
			// const auto& circumfs {
			// 	#ifdef RECTMESH_PRECALC_CIRCUMFS
			// 	MAP( getReadRow( mesh.cells().circumfs() ) )
			// 	#else
			// 	getReadRow( mesh.cells().circumfs() )
			// 	#endif
			// };
			// const auto& areas {
			// 	block.template writeBlock<1>( mesh.cells().areas() )
			// 	.col(meshRow)
			// };
			// const auto& circumfs {
			// 	block.template writeBlock<1>( mesh.cells().circumfs() )
			// 	.col(meshRow)
			// };

			// return 4. * areas / circumfs;

			Eigen::VectorBlock<const ArrayXs> dx {
				mesh.dx().segment( block.xBeg(), block.nx() )
			};
			scalar dy { mesh.dy()[block.yBeg() + y] };
			return 4 * (dx * dy / ( 2 * (dx + dy) )).transpose();
		}
	};

	scalar step { std::numeric_limits<scalar>::max() };

	block.forEachRow(
		[&](Index y){
			step = std::min( step, (
				courantFactor(y) / (
					block.template rowMap<2>( vBlock, y ).colwise().norm() +
					(block.template rowMap<1>( dBlock, y ) * phys::g).sqrt()
				) ).minCoeff()
			);
		}
	);

	return solver().settings().maxCourant * step;
}

} // namespace hms::swe